
##load data
x.vars = c('non.white','hispanic','college.grad','male')
money.dat = read.csv('money_dat.csv')
attribution.dat = read.csv('attribution_dat.csv')
social.dat = read.csv('social_dat.csv')
physical.dat = read.csv('physical_dat_notrim.csv')


cat('\n \n Starting Experiment 2 Estimation, no trimming \n \n') 

##outmatrix
outmat = matrix(nrow = 4, ncol = 11)
rownames(outmat) = c('allocation','attribution','social.perceptions','physical.perceptions')
colnames(outmat) = c('int.mean','int.sd','seg.mean','seg.sd','ate','se','ci.low','ci.high','p.value','cohens.d','N')

###allocation analysis
##RI analysis
money.Xs = as.matrix(money.dat[,x.vars])  ##get necessary x variables
##generate RI permutations and probabilities
money.ri.perms = genperms(money.dat$treatment.assignment,
                        blockvar = money.dat$Block.ID)
money.ri.probs = genprob(money.ri.perms)

##estimate treatment effect and inference
ate = estate(Y = money.dat$allocated.diff,
             Z= money.dat$treatment.assignment,
             prob = money.ri.probs,
             X = money.Xs)
money.Ys = genouts(Y = money.dat$allocated.diff,
                   Z= money.dat$treatment.assignment,
                   ate = 0)
##generate outcomes
distout <- gendist(money.Ys,money.ri.perms, prob=money.ri.probs, X = money.Xs)
money.Ys.ci = genouts(Y = money.dat$allocated.diff,
                      Z= money.dat$treatment.assignment,
                      ate = ate)
distout.ci <- gendist(money.Ys.ci,money.ri.perms, prob=money.ri.probs, X = money.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

##cohen's D for stadndardized effects
cohen = cohen.d(d = money.dat$allocated.diff,
                f = as.factor(money.dat$treatment.assignment))

##add outcomes to matrix
outmat['allocation','int.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 0])
outmat['allocation','int.sd'] = sd(money.dat$allocated.diff[money.dat$treatment.assignment == 0])
outmat['allocation','seg.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 1])
outmat['allocation','seg.sd'] = sd(money.dat$allocated.diff[money.dat$treatment.assignment == 1])
outmat['allocation','ate'] = ate
outmat['allocation','se'] = sd(distout.ci)
outmat['allocation','ci.low'] = quantile(distout.ci,.025)
outmat['allocation','ci.high'] = quantile(distout.ci,.975)
outmat['allocation','p.value'] = disp$greater.p.value
outmat['allocation','cohens.d'] = cohen$estimate
outmat['allocation','N'] = nrow(money.dat)


#################################
##now do same thing for all other outcomes
#################################

#####attribution analysis
attribution.Xs = as.matrix(attribution.dat[,x.vars])
attribution.ri.perms = genperms(attribution.dat$treatment.assignment,
                                blockvar = attribution.dat$Block.ID)
attribution.ri.probs = genprob(attribution.ri.perms)
ate = estate(Y = attribution.dat$attribution.diff,
             Z= attribution.dat$treatment.assignment,
             prob = attribution.ri.probs,
             X = attribution.Xs)
attribution.Ys = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = 0)
distout <- gendist(attribution.Ys,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
attribution.Ys.ci = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = ate)
distout.ci <- gendist(attribution.Ys.ci,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = attribution.dat$attribution.diff,
                f = as.factor(attribution.dat$treatment.assignment))

outmat['attribution','int.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 0])
outmat['attribution','int.sd'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 0])
outmat['attribution','seg.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 1])
outmat['attribution','seg.sd'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 1])
outmat['attribution','ate'] = ate
outmat['attribution','se'] = sd(distout.ci)
outmat['attribution','ci.low'] = quantile(distout.ci,.025)
outmat['attribution','ci.high'] = quantile(distout.ci,.975)
outmat['attribution','p.value'] = disp$greater.p.value
outmat['attribution','cohens.d'] = cohen$estimate
outmat['attribution','N'] = nrow(attribution.dat)


####social analysis
diff.Xs = as.matrix(social.dat[,x.vars])
diff.ri.perms = genperms(social.dat$treatment.assignment,
                         blockvar = social.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = social.dat$social.diff.scale,
             Z= social.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = social.dat$social.diff.scale,
                f = as.factor(social.dat$treatment.assignment))

outmat['social.perceptions','int.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 0])
outmat['social.perceptions','int.sd'] = sd(social.dat$social.diff.scale[social.dat$treatment.assignment == 0])
outmat['social.perceptions','seg.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 1])
outmat['social.perceptions','seg.sd'] = sd(social.dat$social.diff.scale[social.dat$treatment.assignment == 1])
outmat['social.perceptions','ate'] = ate
outmat['social.perceptions','se'] = sd(distout.ci)
outmat['social.perceptions','ci.low'] = quantile(distout.ci,.025)
outmat['social.perceptions','ci.high'] = quantile(distout.ci,.975)
outmat['social.perceptions','p.value'] = disp$greater.p.value
outmat['social.perceptions','cohens.d'] = cohen$estimate
outmat['social.perceptions','N'] = nrow(social.dat)



###physical analysis
diff.Xs = as.matrix(physical.dat[,x.vars])
diff.ri.perms = genperms(physical.dat$treatment.assignment,
                         blockvar = physical.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = physical.dat$physical.diff.scale,
             Z= physical.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = physical.dat$physical.diff.scale,
                f = as.factor(physical.dat$treatment.assignment))


outmat['physical.perceptions','int.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 0])
outmat['physical.perceptions','int.sd'] = sd(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 0])
outmat['physical.perceptions','seg.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 1])
outmat['physical.perceptions','seg.sd'] = sd(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 1])
outmat['physical.perceptions','ate'] = ate
outmat['physical.perceptions','se'] = sd(distout.ci)
outmat['physical.perceptions','ci.low'] = quantile(distout.ci,.025)
outmat['physical.perceptions','ci.high'] = quantile(distout.ci,.975)
outmat['physical.perceptions','p.value'] = disp$greater.p.value
outmat['physical.perceptions','cohens.d'] = cohen$estimate
outmat['physical.perceptions','N'] = nrow(physical.dat)


##format the matrix for output
outmat = as.data.frame(round(outmat,3))
outmat[,1] = paste(outmat[,1], ' (',outmat[,2],')',sep='')
outmat[,3] = paste(outmat[,3], ' (',outmat[,4],')',sep='')
outmat[,5] = paste(outmat[,5], ' (',outmat[,6],')',sep='')
outmat[,7] = paste('[',outmat[,7],',',outmat[,8],']',sep='')

outmat[,8] = NULL; outmat[,6] = NULL; outmat[,4] = NULL; outmat[,2] = NULL

outmat[,'cohens.d'] = outmat[,'cohens.d']*(-1)


##output Appendix Table 2
outmat.reduced = outmat[,c('ate','ci.low','p.value','cohens.d','N')]
table.s12 = xtable(outmat.reduced)
rownames(table.s12) = c('Allocation','Attribution', "Social Perceptions", "Physical Perceptions")
colnames(table.s12) = c('Beta (SE)','CI','P','D','N')

print(table.s12,
      "TableS12.tex",
      type = 'latex',
      floating = F,
      booktabs = T)



cat('\n \n Experiment 2 Results, no trimming: \n')
print(outmat.reduced)



